#!/bin/bash

#=============================================================================
# 基因组组装流水线脚本 (Genome Assembly Pipeline)
# 适用于 PacBio HiFi 长读长测序数据的基因组组装
# 作者: lxh
# 版本: 1.0
# 日期: $(date +%Y-%m-%d)
#=============================================================================

# 设置脚本在遇到错误时退出
set -euo pipefail

# 颜色输出函数
RED='\033[0;31m'
GREEN='\033[0;32m'
YELLOW='\033[1;33m'
BLUE='\033[0;34m'
NC='\033[0m' # No Color

log_info() { echo -e "${GREEN}[INFO]${NC} $1"; }
log_warn() { echo -e "${YELLOW}[WARN]${NC} $1"; }
log_error() { echo -e "${RED}[ERROR]${NC} $1"; }
log_step() { echo -e "${BLUE}[STEP]${NC} $1"; }

# 显示帮助信息
show_help() {
    cat << EOF
基因组组装流水线脚本

用法: $0 [选项]

必需参数:
  -i, --input FILE        输入的BAM文件路径
  -p, --prefix STRING     输出文件前缀

可选参数:
  -t, --threads INT       线程数 (默认: 8)
  -m, --memory INT        内存大小(GB) (默认: 48)
  -k, --kmer INT          k-mer大小 (默认: 21)
  -n, --ploidy INT        倍性 (1=单倍体, 2=二倍体) (默认: 1)
  --telo STRING           端粒序列 (默认: TAACCC)
  --hifiasm-k INT         hifiasm k-mer大小 (默认: 63)
  --hifiasm-r INT         hifiasm rounds数 (默认: 3)
  --hifiasm-l INT         hifiasm 长度阈值 (默认: 0)
  --samtools-threads INT  samtools线程数 (默认: 使用全局线程数)
  --kmc-threads INT       kmc线程数 (默认: 使用全局线程数)
  --hifiasm-threads INT   hifiasm线程数 (默认: 使用全局线程数)
  --kmc-ci INT            kmc最小计数阈值 (默认: 3)
  --kmc-cs INT            kmc最大计数阈值 (默认: 10000)
  --skip-md5              跳过MD5验证
  --skip-survey           跳过基因组调查
  --skip-assembly         跳过基因组组装
  -h, --help              显示此帮助信息

示例:
  $0 -i sample.bam -p my_genome -t 16 -m 64 -n 2
  $0 --input data.bam --prefix genome_v1 --threads 12 --memory 32 --ploidy 1
  $0 -i hifi.bam -p assembly --hifiasm-k 51 --hifiasm-r 5 --hifiasm-threads 16
  $0 -i data.bam -p test --kmc-ci 5 --kmc-cs 50000 --kmc-threads 8

EOF
}

# 默认参数设置
INPUT_BAM=""
PREFIX=""
THREADS=8
MEMORY=48
KMER_SIZE=21
PLOIDY=1
TELO_SEQ="TAACCC"
HIFIASM_K=63
HIFIASM_R=3
HIFIASM_L=0
SAMTOOLS_THREADS=""
KMC_THREADS=""
HIFIASM_THREADS=""
KMC_CI=3
KMC_CS=10000
SKIP_MD5=false
SKIP_SURVEY=false
SKIP_ASSEMBLY=false

# 解析命令行参数
while [[ $# -gt 0 ]]; do
    case $1 in
        -i|--input)
            INPUT_BAM="$2"
            shift 2
            ;;
        -p|--prefix)
            PREFIX="$2"
            shift 2
            ;;
        -t|--threads)
            THREADS="$2"
            shift 2
            ;;
        -m|--memory)
            MEMORY="$2"
            shift 2
            ;;
        -k|--kmer)
            KMER_SIZE="$2"
            shift 2
            ;;
        -n|--ploidy)
            PLOIDY="$2"
            shift 2
            ;;
        --telo)
            TELO_SEQ="$2"
            shift 2
            ;;
        --hifiasm-k)
            HIFIASM_K="$2"
            shift 2
            ;;
        --hifiasm-r)
            HIFIASM_R="$2"
            shift 2
            ;;
        --hifiasm-l)
            HIFIASM_L="$2"
            shift 2
            ;;
        --samtools-threads)
            SAMTOOLS_THREADS="$2"
            shift 2
            ;;
        --kmc-threads)
            KMC_THREADS="$2"
            shift 2
            ;;
        --hifiasm-threads)
            HIFIASM_THREADS="$2"
            shift 2
            ;;
        --kmc-ci)
            KMC_CI="$2"
            shift 2
            ;;
        --kmc-cs)
            KMC_CS="$2"
            shift 2
            ;;
        --skip-md5)
            SKIP_MD5=true
            shift
            ;;
        --skip-survey)
            SKIP_SURVEY=true
            shift
            ;;
        --skip-assembly)
            SKIP_ASSEMBLY=true
            shift
            ;;
        -h|--help)
            show_help
            exit 0
            ;;
        *)
            log_error "未知参数: $1"
            show_help
            exit 1
            ;;
    esac
done

# 验证必需参数
if [[ -z "$INPUT_BAM" || -z "$PREFIX" ]]; then
    log_error "缺少必需参数！"
    show_help
    exit 1
fi

# 验证输入文件存在
if [[ ! -f "$INPUT_BAM" ]]; then
    log_error "输入文件不存在: $INPUT_BAM"
    exit 1
fi

# 验证参数范围
if [[ $THREADS -lt 1 || $THREADS -gt 128 ]]; then
    log_error "线程数必须在 1-128 之间"
    exit 1
fi

if [[ $MEMORY -lt 1 || $MEMORY -gt 1024 ]]; then
    log_error "内存大小必须在 1-1024 GB 之间"
    exit 1
fi

if [[ $PLOIDY -ne 1 && $PLOIDY -ne 2 ]]; then
    log_error "倍性只能是 1 (单倍体) 或 2 (二倍体)"
    exit 1
fi

if [[ $HIFIASM_K -lt 1 || $HIFIASM_K -gt 127 ]]; then
    log_error "hifiasm k-mer大小必须在 1-127 之间"
    exit 1
fi

if [[ $HIFIASM_R -lt 1 || $HIFIASM_R -gt 10 ]]; then
    log_error "hifiasm rounds数必须在 1-10 之间"
    exit 1
fi

if [[ $HIFIASM_L -lt 0 ]]; then
    log_error "hifiasm 长度阈值不能为负数"
    exit 1
fi

if [[ $KMC_CI -lt 1 || $KMC_CI -gt 1000 ]]; then
    log_error "kmc最小计数阈值必须在 1-1000 之间"
    exit 1
fi

if [[ $KMC_CS -lt 100 || $KMC_CS -gt 100000 ]]; then
    log_error "kmc最大计数阈值必须在 100-100000 之间"
    exit 1
fi

# 设置默认线程数（如果未指定）
if [[ -z "$SAMTOOLS_THREADS" ]]; then
    SAMTOOLS_THREADS="$THREADS"
fi

if [[ -z "$KMC_THREADS" ]]; then
    KMC_THREADS="$THREADS"
fi

if [[ -z "$HIFIASM_THREADS" ]]; then
    HIFIASM_THREADS="$THREADS"
fi

# 检查必需的软件工具并记录版本信息
check_tools() {
    local tools=("samtools" "seqkit" "kmc" "kmc_tools" "hifiasm")
    local missing_tools=()
    
    # 常见的工具安装路径
    local common_paths=(
        "/home/$USER/Opt/compile/Samtools_v1.22.1"
        "/home/$USER/Opt/bin"
        "/home/$USER/Opt/src/kmc_v3.2.4/bin"
        "/home/$USER/Opt/compile/Hifiasm_v0.25.0"
        "/usr/local/bin"
        "/opt/bin"
        "$HOME/bin"
        "$HOME/.local/bin"
    )
    
    # 添加常见路径到PATH (如果不存在)
    for path in "${common_paths[@]}"; do
        if [[ -d "$path" && ":$PATH:" != *":$path:"* ]]; then
            export PATH="$path:$PATH"
        fi
    done
    
    # 创建临时目录（如果不存在）
    mkdir -p "${PREFIX}_temp"
    
    # 创建版本信息文件
    local version_file="${PREFIX}_temp/software_versions.txt"
    echo "软件版本信息" > "$version_file"
    echo "==================" >> "$version_file"
    echo "检查时间: $(date)" >> "$version_file"
    echo "" >> "$version_file"
    
    for tool in "${tools[@]}"; do
        if ! command -v "$tool" &> /dev/null; then
            missing_tools+=("$tool")
        else
            local tool_path=$(command -v "$tool")
            log_info "找到工具: $tool -> $tool_path"
            
            # 获取软件版本信息
            case $tool in
                "samtools")
                    local version=$(samtools --version 2>/dev/null | head -1 | cut -d' ' -f2 2>/dev/null)
                    if [[ -z "$version" || "$version" == "未知版本" ]]; then
                        # 尝试从路径提取版本信息
                        version=$(echo "$tool_path" | grep -o "v[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1 | sed 's/^v//')
                        [[ -z "$version" ]] && version=$(echo "$tool_path" | grep -o "[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1)
                        [[ -z "$version" ]] && version="路径未包含版本信息"
                    fi
                    echo "$tool: $version ($tool_path)" >> "$version_file"
                    ;;
                "seqkit")
                    local version=$(seqkit version 2>/dev/null | grep "seqkit" | cut -d' ' -f2 2>/dev/null | sed 's/^v//')
                    if [[ -z "$version" || "$version" == "未知版本" ]]; then
                        version=$(echo "$tool_path" | grep -o "v[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1 | sed 's/^v//')
                        [[ -z "$version" ]] && version=$(echo "$tool_path" | grep -o "[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1)
                        [[ -z "$version" ]] && version="路径未包含版本信息"
                    fi
                    echo "$tool: $version ($tool_path)" >> "$version_file"
                    ;;
                "kmc")
                    local version=$(kmc 2>&1 | head -1 | grep -o "v[0-9]\.[0-9]\.[0-9]" 2>/dev/null | sed 's/^v//')
                    if [[ -z "$version" || "$version" == "未知版本" ]]; then
                        version=$(echo "$tool_path" | grep -o "v[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1 | sed 's/^v//')
                        [[ -z "$version" ]] && version=$(echo "$tool_path" | grep -o "[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1)
                        [[ -z "$version" ]] && version="3.2.4(从路径推断)"
                    fi
                    echo "$tool: $version ($tool_path)" >> "$version_file"
                    ;;
                "kmc_tools")
                    local version=$(kmc_tools 2>&1 | head -1 | grep -o "v[0-9]\.[0-9]\.[0-9]" 2>/dev/null | sed 's/^v//')
                    if [[ -z "$version" || "$version" == "未知版本" ]]; then
                        version=$(echo "$tool_path" | grep -o "v[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1 | sed 's/^v//')
                        [[ -z "$version" ]] && version=$(echo "$tool_path" | grep -o "[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1)
                        [[ -z "$version" ]] && version="3.2.4(从路径推断)"
                    fi
                    echo "$tool: $version ($tool_path)" >> "$version_file"
                    ;;
                "hifiasm")
                    local version=$(hifiasm --version 2>&1 | head -1 | cut -d'-' -f1 2>/dev/null)
                    if [[ -z "$version" || "$version" == "未知版本" ]]; then
                        version=$(echo "$tool_path" | grep -o "v[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1 | sed 's/^v//')
                        [[ -z "$version" ]] && version=$(echo "$tool_path" | grep -o "[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1)
                        [[ -z "$version" ]] && version="0.25.0(从路径推断)"
                    fi
                    echo "$tool: $version ($tool_path)" >> "$version_file"
                    ;;
            esac
        fi
    done
    
    if [[ ${#missing_tools[@]} -gt 0 ]]; then
        log_error "缺少以下工具: ${missing_tools[*]}"
        log_error "当前PATH: $PATH"
        log_error "请确保这些工具在PATH中或使用 'sudo -E' 运行脚本"
        exit 1
    fi
    
    # 检查可选工具
    if ! command -v "genomescope.R" &> /dev/null && [[ $SKIP_SURVEY == false ]]; then
        log_warn "genomescope.R 未找到，将跳过基因组调查步骤"
        SKIP_SURVEY=true
    else
        local genomescope_path=$(command -v genomescope.R)
        local genomescope_version=$(Rscript -e "packageVersion('GenomeScope')" 2>/dev/null)
        
        # 如果无法从R包获取版本，尝试从路径提取
        if [[ -z "$genomescope_version" || "$genomescope_version" == "未知版本" ]]; then
            genomescope_version=$(echo "$genomescope_path" | grep -o "v[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1 | sed 's/^v//')
            [[ -z "$genomescope_version" ]] && genomescope_version=$(echo "$genomescope_path" | grep -o "[0-9]\+\.[0-9]\+\.[0-9]\+" | head -1)
            [[ -z "$genomescope_version" ]] && genomescope_version="2.0.1(从路径推断)"
        fi
        
        echo "genomescope.R: $genomescope_version ($genomescope_path)" >> "$version_file"
    fi
    
    log_info "软件版本信息已保存到: $version_file"
}

# 创建输出目录
create_dirs() {
    local dirs=("${PREFIX}_logs" "${PREFIX}_temp" "${PREFIX}_results")
    for dir in "${dirs[@]}"; do
        mkdir -p "$dir"
        log_info "创建目录: $dir"
    done
    
    # 创建子目录用于组织文件
    mkdir -p "${PREFIX}_temp/kmc_files"
    mkdir -p "${PREFIX}_temp/assembly_files"
    mkdir -p "${PREFIX}_results/genomescope"
    mkdir -p "${PREFIX}_results/stats"
}

# 记录运行参数
log_parameters() {
    local log_file="${PREFIX}_logs/parameters.log"
    cat > "$log_file" << EOF
基因组组装流水线运行参数
========================
运行时间: $(date)
输入文件: $INPUT_BAM
输出前缀: $PREFIX
线程数: $THREADS
内存大小: ${MEMORY}GB
K-mer大小: $KMER_SIZE
倍性: $PLOIDY
端粒序列: $TELO_SEQ
Hifiasm K-mer大小: $HIFIASM_K
Hifiasm Rounds数: $HIFIASM_R
Hifiasm 长度阈值: $HIFIASM_L
Samtools线程数: $SAMTOOLS_THREADS
KMC线程数: $KMC_THREADS
Hifiasm线程数: $HIFIASM_THREADS
KMC最小计数阈值: $KMC_CI
KMC最大计数阈值: $KMC_CS
跳过MD5验证: $SKIP_MD5
跳过基因组调查: $SKIP_SURVEY
跳过基因组组装: $SKIP_ASSEMBLY
========================
EOF
    log_info "运行参数已保存到: $log_file"
}

# 步骤1: 数据提取和统计
step1_data_extraction() {
    log_step "步骤1: 数据提取和统计"
    
    # MD5验证 (如果存在md5sum.txt文件且未跳过)
    if [[ -f "md5sum.txt" && $SKIP_MD5 == false ]]; then
        log_info "验证文件完整性..."
        if md5sum -c md5sum.txt; then
            log_info "MD5验证通过"
        else
            log_error "MD5验证失败，请检查文件完整性"
            exit 1
        fi
    else
        log_warn "跳过MD5验证"
    fi
    
    # 从BAM文件提取FASTQ数据
    local output_fastq="${PREFIX}_temp/${PREFIX}_hifiread_all.fastq.gz"
    log_info "从BAM文件提取FASTQ数据..."
    log_info "输入: $INPUT_BAM"
    log_info "输出: $output_fastq"
    
    samtools fastq -@ "$SAMTOOLS_THREADS" "$INPUT_BAM" | \
        seqkit seq -j "$SAMTOOLS_THREADS" -w 0 -o "$output_fastq"
    
    if [[ ! -f "$output_fastq" ]]; then
        log_error "FASTQ文件生成失败"
        exit 1
    fi
    
    # 生成reads统计信息
    log_info "生成reads统计信息..."
    seqkit stats -j "$SAMTOOLS_THREADS" -a -o "${PREFIX}_results/stats/${PREFIX}_hifiread_stats.xls" "$output_fastq"
    
    # 生成reads长度分布直方图
    log_info "生成reads长度分布直方图..."
    seqkit watch -j "$SAMTOOLS_THREADS" -f ReadLen -O "${PREFIX}_results/stats/${PREFIX}_hifiread_length_histogram.png" "$output_fastq"
    
    log_info "步骤1完成"
}

# 步骤2: 基因组调查
step2_genome_survey() {
    if [[ $SKIP_SURVEY == true ]]; then
        log_warn "跳过基因组调查步骤"
        return
    fi
    
    log_step "步骤2: 基因组调查"
    
    local input_fastq="${PREFIX}_temp/${PREFIX}_hifiread_all.fastq.gz"
    local kmc_output="${PREFIX}_temp/kmc_files/${PREFIX}_kmc"
    local histo_file="${PREFIX}_temp/kmc_files/${PREFIX}_${KMER_SIZE}kmer.histo"
    
    # k-mer计数分析
    log_info "进行k-mer计数分析 (k=$KMER_SIZE)..."
    kmc -k"$KMER_SIZE" -m"$MEMORY" -t"$KMC_THREADS" -ci"$KMC_CI" -cs"$KMC_CS" \
        "$input_fastq" "$kmc_output" "${PREFIX}_temp/kmc_files/"
    
    # 生成k-mer频率直方图
    log_info "生成k-mer频率直方图..."
    kmc_tools transform "$kmc_output" histogram "$histo_file" -cx10000
    
    # 基因组特征估算 (如果genomescope.R可用)
    if command -v "genomescope.R" &> /dev/null; then
        log_info "使用GenomeScope估算基因组特征..."
        genomescope.R -i "$histo_file" -o "${PREFIX}_results/genomescope" -k "$KMER_SIZE" -p "$PLOIDY"
    else
        log_warn "genomescope.R不可用，跳过基因组特征估算"
    fi
    
    log_info "步骤2完成"
}

# 步骤3: 基因组组装
step3_genome_assembly() {
    if [[ $SKIP_ASSEMBLY == true ]]; then
        log_warn "跳过基因组组装步骤"
        return
    fi
    
    log_step "步骤3: 基因组组装"
    
    local input_fastq="${PREFIX}_temp/${PREFIX}_hifiread_all.fastq.gz"
    local assembly_prefix="${PREFIX}_temp/assembly_files/${PREFIX}_hifiasm_run001"
    
    # 设置hifiasm参数
    local hifiasm_params="-o $assembly_prefix -k $HIFIASM_K -r $HIFIASM_R -l $HIFIASM_L --n-hap $PLOIDY -t $HIFIASM_THREADS"
    if [[ -n "$TELO_SEQ" ]]; then
        hifiasm_params="$hifiasm_params --telo-m $TELO_SEQ"
    fi
    
    # 运行hifiasm组装
    log_info "使用hifiasm进行基因组组装..."
    log_info "参数: $hifiasm_params"
    hifiasm $hifiasm_params "$input_fastq"
    
    # 转换GFA为FASTA格式
    log_info "转换GFA格式为FASTA格式..."
    local gfa_file="${assembly_prefix}.bp.p_ctg.gfa"
    local fasta_file="${PREFIX}_results/${PREFIX}_PCTG.fasta"
    local compressed_fasta="${fasta_file}.gz"
    
    if [[ -f "$gfa_file" ]]; then
        # 提取contigs序列
        awk -F "[\t:]" '{if($1=="S"){printf ">%s\tlength=%s\treads=%s\n%s\n",$2,$6,$9,$3;}}' \
            "$gfa_file" | tee "$fasta_file" | seqkit seq -w 100 -o "$compressed_fasta"
        
        # 生成组装统计信息
        log_info "生成组装统计信息..."
        seqkit stats -j "$SAMTOOLS_THREADS" -a -o "${PREFIX}_results/stats/${PREFIX}_PCTG.fasta.stats.xls" "$compressed_fasta"
        
        # 提取contigs连接信息
        log_info "提取contigs连接信息..."
        awk '{if($1=="L"){print;}}' "$gfa_file" > "${PREFIX}_results/${PREFIX}_PCTG.links"
        
        log_info "组装完成，结果文件保存在 ${PREFIX}_results/ 目录"
    else
        log_error "未找到GFA文件: $gfa_file"
        exit 1
    fi
    
    log_info "步骤3完成"
}

# 生成markdown格式运行报告
generate_report() {
    log_step "生成运行报告"
    
    local report_file="${PREFIX}_results/${PREFIX}_assembly_report.md"
    local start_time=$(cat "${PREFIX}_logs/parameters.log" | grep "运行时间" | cut -d' ' -f2-)
    local end_time=$(date)
    
    cat > "$report_file" << 'EOF'
# 基因组组装流水线运行报告

## 基本信息

**运行时间**: %START_TIME%  
**完成时间**: %END_TIME%  
**输出前缀**: %PREFIX%  
**输入文件**: %INPUT_BAM%  

## 运行参数

| 参数 | 值 |
|------|-----|
| 线程数 | %THREADS% |
| 内存大小 | %MEMORY%GB |
| K-mer大小 | %KMER_SIZE% |
| 倍性 | %PLOIDY% |
| 端粒序列 | %TELO_SEQ% |
| Hifiasm K-mer大小 | %HIFIASM_K% |
| Hifiasm Rounds数 | %HIFIASM_R% |
| Hifiasm 长度阈值 | %HIFIASM_L% |
| Samtools线程数 | %SAMTOOLS_THREADS% |
| KMC线程数 | %KMC_THREADS% |
| Hifiasm线程数 | %HIFIASM_THREADS% |
| KMC最小计数阈值 | %KMC_CI% |
| KMC最大计数阈值 | %KMC_CS% |

## 软件版本信息

```
%SOFTWARE_VERSIONS%
```

## 数据统计结果

### 1. HiFi Reads统计

EOF
    
    # 插入实际值
    sed -i "s|%START_TIME%|$start_time|g" "$report_file"
    sed -i "s|%END_TIME%|$end_time|g" "$report_file"
    sed -i "s|%PREFIX%|$PREFIX|g" "$report_file"
    sed -i "s|%INPUT_BAM%|$INPUT_BAM|g" "$report_file"
    sed -i "s|%THREADS%|$THREADS|g" "$report_file"
    sed -i "s|%MEMORY%|$MEMORY|g" "$report_file"
    sed -i "s|%KMER_SIZE%|$KMER_SIZE|g" "$report_file"
    sed -i "s|%PLOIDY%|$PLOIDY|g" "$report_file"
    sed -i "s|%TELO_SEQ%|$TELO_SEQ|g" "$report_file"
    sed -i "s|%HIFIASM_K%|$HIFIASM_K|g" "$report_file"
    sed -i "s|%HIFIASM_R%|$HIFIASM_R|g" "$report_file"
    sed -i "s|%HIFIASM_L%|$HIFIASM_L|g" "$report_file"
    sed -i "s|%SAMTOOLS_THREADS%|$SAMTOOLS_THREADS|g" "$report_file"
    sed -i "s|%KMC_THREADS%|$KMC_THREADS|g" "$report_file"
    sed -i "s|%HIFIASM_THREADS%|$HIFIASM_THREADS|g" "$report_file"
    sed -i "s|%KMC_CI%|$KMC_CI|g" "$report_file"
    sed -i "s|%KMC_CS%|$KMC_CS|g" "$report_file"
    
    # 插入软件版本信息
    # 首先尝试从results目录读取（如果已经备份），否则从temp目录读取
    local version_file=""
    if [[ -f "${PREFIX}_results/software_versions.txt" ]]; then
        version_file="${PREFIX}_results/software_versions.txt"
    elif [[ -f "${PREFIX}_temp/software_versions.txt" ]]; then
        version_file="${PREFIX}_temp/software_versions.txt"
    fi
    
    if [[ -n "$version_file" && -f "$version_file" ]]; then
        # 创建临时文件来逐步构建最终报告
        local temp_report=$(mktemp)
        
        # 读取软件版本信息到变量，转义特殊字符
        local versions_content=""
        while IFS= read -r line; do
            versions_content="${versions_content}${line}\n"
        done < "$version_file"
        
        # 使用awk进行替换，避免sed的多行问题
        awk -v versions="$versions_content" '
        /%SOFTWARE_VERSIONS%/ {
            printf "%s", versions
            next
        }
        { print }
        ' "$report_file" > "$temp_report"
        
        mv "$temp_report" "$report_file"
    fi
    
    # 添加HiFi reads统计表格
    if [[ -f "${PREFIX}_results/stats/${PREFIX}_hifiread_stats.xls" ]]; then
        echo "" >> "$report_file"
        
        # 检查是否能正确解析统计文件
        if [[ $(wc -l < "${PREFIX}_results/stats/${PREFIX}_hifiread_stats.xls") -gt 1 ]]; then
            # 尝试解析数据
            local stats_line=$(tail -n +2 "${PREFIX}_results/stats/${PREFIX}_hifiread_stats.xls" | head -1)
            if [[ -n "$stats_line" ]]; then
                # 使用数组读取数据
                IFS=$'\t' read -ra STATS_ARRAY <<< "$stats_line"
                
                # 检查数组长度是否足够
                if [[ ${#STATS_ARRAY[@]} -ge 17 ]]; then
                    echo "| 指标 | 值 |" >> "$report_file"
                    echo "|------|-----|" >> "$report_file"
                    echo "| 序列数量 | ${STATS_ARRAY[3]} |" >> "$report_file"
                    echo "| 总长度 | ${STATS_ARRAY[4]} bp |" >> "$report_file"
                    echo "| 平均长度 | ${STATS_ARRAY[6]} bp |" >> "$report_file"
                    echo "| 最小长度 | ${STATS_ARRAY[5]} bp |" >> "$report_file"
                    echo "| 最大长度 | ${STATS_ARRAY[7]} bp |" >> "$report_file"
                    echo "| N50 | ${STATS_ARRAY[11]} bp |" >> "$report_file"
                    echo "| Q20质量 | ${STATS_ARRAY[13]}% |" >> "$report_file"
                    echo "| Q30质量 | ${STATS_ARRAY[14]}% |" >> "$report_file"
                    echo "| 平均质量 | ${STATS_ARRAY[15]} |" >> "$report_file"
                    echo "| GC含量 | ${STATS_ARRAY[16]}% |" >> "$report_file"
                else
                    # 如果解析失败，且原始文件非空，则仅输出原始内容（不保留空标题）
                    if [[ -s "${PREFIX}_results/stats/${PREFIX}_hifiread_stats.xls" ]]; then
                        echo "```" >> "$report_file"
                        cat "${PREFIX}_results/stats/${PREFIX}_hifiread_stats.xls" >> "$report_file"
                        echo "```" >> "$report_file"
                    else
                        log_warn "HiFi reads统计原始文件无内容，跳过原始结果显示"
                    fi
                fi
            else
                # 无法读取数据时不添加空标题，跳过整个部分
                log_warn "无法读取HiFi reads统计数据，跳过统计表格显示"
            fi
        else
            # 文件为空时不添加空标题，跳过整个部分  
            log_warn "HiFi reads统计文件为空，跳过统计表格显示"
        fi
        
        # 添加reads长度分布直方图
        echo "" >> "$report_file"
        echo "#### HiFi Reads长度分布" >> "$report_file"
        echo "" >> "$report_file"
        if [[ -f "${PREFIX}_results/stats/${PREFIX}_hifiread_length_histogram.png" ]]; then
            echo "![HiFi Reads长度分布直方图](stats/${PREFIX}_hifiread_length_histogram.png)" >> "$report_file"
        else
            echo "*注：长度分布直方图文件未找到*" >> "$report_file"
        fi
    else
        echo "" >> "$report_file"
        echo "*HiFi reads统计文件未找到*" >> "$report_file"
    fi
    
    # 添加k-mer统计结果
    echo "" >> "$report_file"
    echo "### 2. K-mer分析结果" >> "$report_file"
    echo "" >> "$report_file"
    
    # 添加k-mer频率直方图信息
    local kmer_histo_file=""
    # 首先检查temp目录，如果没有则检查results目录中是否有备份
    if [[ -f "${PREFIX}_temp/kmc_files/${PREFIX}_${KMER_SIZE}kmer.histo" ]]; then
        kmer_histo_file="${PREFIX}_temp/kmc_files/${PREFIX}_${KMER_SIZE}kmer.histo"
    elif [[ -f "${PREFIX}_results/kmc_files/${PREFIX}_${KMER_SIZE}kmer.histo" ]]; then
        kmer_histo_file="${PREFIX}_results/kmc_files/${PREFIX}_${KMER_SIZE}kmer.histo"
    fi
    
    if [[ -n "$kmer_histo_file" && -f "$kmer_histo_file" ]]; then
        echo "#### K-mer频率分布" >> "$report_file"
        echo "" >> "$report_file"
        echo "- **K-mer大小**: ${KMER_SIZE}" >> "$report_file"
        echo "- **频率直方图文件**: $(basename "$kmer_histo_file")" >> "$report_file"
        
        # 统计k-mer数量
        local total_kmers=$(awk '{sum+=$2} END {print sum}' "$kmer_histo_file" 2>/dev/null || echo "未知")
        local unique_kmers=$(wc -l < "$kmer_histo_file" 2>/dev/null || echo "未知")
        
        echo "" >> "$report_file"
        echo "| K-mer统计指标 | 值 |" >> "$report_file"
        echo "|--------------|-----|" >> "$report_file"
        echo "| 总K-mer数量 | $total_kmers |" >> "$report_file"
        echo "| 唯一K-mer数量 | $unique_kmers |" >> "$report_file"
    else
        echo "#### K-mer频率分布" >> "$report_file"
        echo "" >> "$report_file"
        echo "- **K-mer大小**: ${KMER_SIZE}" >> "$report_file"
        echo "- **状态**: K-mer直方图文件已清理或未生成" >> "$report_file"
    fi
    
    # 添加基因组调查结果（summary与图表在未跳过时显示；model.txt只要存在就显示）
    if [[ ! $SKIP_SURVEY == true ]]; then
        echo "" >> "$report_file"
        echo "### 3. 基因组调查结果 (GenomeScope)" >> "$report_file"
        echo "" >> "$report_file"
        
        # 显示summary.txt内容
        if [[ -f "${PREFIX}_results/genomescope/summary.txt" && -s "${PREFIX}_results/genomescope/summary.txt" ]]; then
            echo "#### 基因组特征估算摘要" >> "$report_file"
            echo "" >> "$report_file"
            
            # 解析summary.txt中的关键信息并格式化显示
            echo "| 基因组特征 | 估算值 |" >> "$report_file"
            echo "|-----------|--------|" >> "$report_file"
            
            # 提取关键数据
            if grep -q "Genome Haploid Length" "${PREFIX}_results/genomescope/summary.txt"; then
                local haploid_length=$(grep "Genome Haploid Length" "${PREFIX}_results/genomescope/summary.txt" | awk '{print $4" bp - "$6" bp"}')
                echo "| 单倍体基因组长度 | $haploid_length |" >> "$report_file"
            fi
            
            if grep -q "Genome Repeat Length" "${PREFIX}_results/genomescope/summary.txt"; then
                local repeat_length=$(grep "Genome Repeat Length" "${PREFIX}_results/genomescope/summary.txt" | awk '{print $4" bp - "$6" bp"}')
                echo "| 重复序列长度 | $repeat_length |" >> "$report_file"
            fi
            
            if grep -q "Genome Unique Length" "${PREFIX}_results/genomescope/summary.txt"; then
                local unique_length=$(grep "Genome Unique Length" "${PREFIX}_results/genomescope/summary.txt" | awk '{print $4" bp - "$6" bp"}')
                echo "| 唯一序列长度 | $unique_length |" >> "$report_file"
            fi
            
            if grep -q "Model Fit" "${PREFIX}_results/genomescope/summary.txt"; then
                local model_fit=$(grep "Model Fit" "${PREFIX}_results/genomescope/summary.txt" | awk '{print $3" - "$4}')
                echo "| 模型拟合度 | $model_fit |" >> "$report_file"
            fi
            
            if grep -q "Read Error Rate" "${PREFIX}_results/genomescope/summary.txt"; then
                local error_rate=$(grep "Read Error Rate" "${PREFIX}_results/genomescope/summary.txt" | awk '{print $4}')
                echo "| 测序错误率 | $error_rate |" >> "$report_file"
            fi
            
            echo "" >> "$report_file"
            echo "<details>" >> "$report_file"
            echo "<summary>完整的GenomeScope摘要信息</summary>" >> "$report_file"
            echo "" >> "$report_file"
            echo "```" >> "$report_file"
            cat "${PREFIX}_results/genomescope/summary.txt" >> "$report_file"
            echo "```" >> "$report_file"
            echo "</details>" >> "$report_file"
            echo "" >> "$report_file"
        fi
        
        # 显示GenomeScope分析图表
        if [[ -f "${PREFIX}_results/genomescope/linear_plot.png" || -f "${PREFIX}_results/genomescope/log_plot.png" ]]; then
            echo "#### GenomeScope分析图表" >> "$report_file"
            echo "" >> "$report_file"
        fi
        
        if [[ -f "${PREFIX}_results/genomescope/linear_plot.png" ]]; then
            echo "**线性图：**" >> "$report_file"
            echo "" >> "$report_file"
            echo "![GenomeScope线性图](genomescope/linear_plot.png)" >> "$report_file"
            echo "" >> "$report_file"
        fi
        
        if [[ -f "${PREFIX}_results/genomescope/log_plot.png" ]]; then
            echo "**对数图：**" >> "$report_file"
            echo "" >> "$report_file"
            echo "![GenomeScope对数图](genomescope/log_plot.png)" >> "$report_file"
            echo "" >> "$report_file"
        fi
        
    fi
    # model.txt 只要存在就显示（即使跳过了GenomeScope计算）
    if [[ -f "${PREFIX}_results/genomescope/model.txt" && -s "${PREFIX}_results/genomescope/model.txt" ]]; then
        echo "" >> "$report_file"
        echo "### GenomeScope模型参数" >> "$report_file"
        echo "" >> "$report_file"
        echo "<details>" >> "$report_file"
        echo "<summary>点击展开查看model.txt</summary>" >> "$report_file"
        echo "" >> "$report_file"
        echo "```" >> "$report_file"
        cat "${PREFIX}_results/genomescope/model.txt" >> "$report_file"
        echo "```" >> "$report_file"
        echo "</details>" >> "$report_file"
        echo "" >> "$report_file"
    fi
    
    # 添加组装结果
    if [[ ! $SKIP_ASSEMBLY == true ]]; then
        echo "" >> "$report_file"
        echo "### 4. 基因组组装结果 (Hifiasm)" >> "$report_file"
        echo "" >> "$report_file"
        
        # 首先显示hifiasm参数信息
        echo "#### Hifiasm组装参数" >> "$report_file"
        echo "" >> "$report_file"
        echo "| 参数 | 值 |" >> "$report_file"
        echo "|------|-----|" >> "$report_file"
        echo "| K-mer大小 | $HIFIASM_K |" >> "$report_file"
        echo "| Rounds数 | $HIFIASM_R |" >> "$report_file"
        echo "| 长度阈值 | $HIFIASM_L |" >> "$report_file"
        echo "| 倍性 | $PLOIDY |" >> "$report_file"
        echo "| 端粒序列 | $TELO_SEQ |" >> "$report_file"
        echo "| 使用线程数 | $HIFIASM_THREADS |" >> "$report_file"
        echo "" >> "$report_file"
        
        # 组装统计结果
        if [[ -f "${PREFIX}_results/stats/${PREFIX}_PCTG.fasta.stats.xls" ]]; then
            echo "#### 主要组装统计 (Primary Contigs)" >> "$report_file"
            echo "" >> "$report_file"
            
            # 检查是否能正确解析组装统计文件
            if [[ $(wc -l < "${PREFIX}_results/stats/${PREFIX}_PCTG.fasta.stats.xls") -gt 1 ]]; then
                # 尝试解析数据
                local assembly_stats_line=$(tail -n +2 "${PREFIX}_results/stats/${PREFIX}_PCTG.fasta.stats.xls" | head -1)
                if [[ -n "$assembly_stats_line" ]]; then
                    # 使用数组读取数据
                    IFS=$'\t' read -ra ASSEMBLY_ARRAY <<< "$assembly_stats_line"
                    
                    # 检查数组长度是否足够
                    if [[ ${#ASSEMBLY_ARRAY[@]} -ge 17 ]]; then
                        echo "| 指标 | 值 |" >> "$report_file"
                        echo "|------|-----|" >> "$report_file"
                        echo "| Contigs数量 | ${ASSEMBLY_ARRAY[3]} |" >> "$report_file"
                        echo "| 总组装长度 | ${ASSEMBLY_ARRAY[4]} bp |" >> "$report_file"
                        echo "| 平均Contig长度 | ${ASSEMBLY_ARRAY[6]} bp |" >> "$report_file"
                        echo "| 最小Contig长度 | ${ASSEMBLY_ARRAY[5]} bp |" >> "$report_file"
                        echo "| 最大Contig长度 | ${ASSEMBLY_ARRAY[7]} bp |" >> "$report_file"
                        echo "| N50长度 | ${ASSEMBLY_ARRAY[11]} bp |" >> "$report_file"
                        echo "| N50对应的contig编号 | ${ASSEMBLY_ARRAY[12]} |" >> "$report_file"
                        echo "| GC含量 | ${ASSEMBLY_ARRAY[16]}% |" >> "$report_file"
                    else
                        # 如果解析失败，且原始文件非空，则仅输出原始内容（不保留空标题）
                        if [[ -s "${PREFIX}_results/stats/${PREFIX}_PCTG.fasta.stats.xls" ]]; then
                            echo "```" >> "$report_file"
                            cat "${PREFIX}_results/stats/${PREFIX}_PCTG.fasta.stats.xls" >> "$report_file"
                            echo "```" >> "$report_file"
                        else
                            log_warn "组装统计原始文件无内容，跳过原始结果显示"
                        fi
                    fi
                else
                    # 无法读取数据时跳过整个部分
                    log_warn "无法读取组装统计数据，跳过统计表格显示"
                fi
            else
                # 文件为空时跳过整个部分
                log_warn "组装统计文件为空，跳过统计表格显示"
            fi
        else
            # 文件不存在时跳过整个部分
            log_warn "组装统计文件未找到，跳过统计表格显示"
        fi
        
        echo "" >> "$report_file"
        echo "#### Hifiasm输出文件说明" >> "$report_file"
        echo "" >> "$report_file"
        echo "| 文件类型 | 文件名 | 说明 |" >> "$report_file"
        echo "|---------|-------|------|" >> "$report_file"
        echo "| Primary contigs (GFA) | ${PREFIX}_hifiasm_run001.bp.p_ctg.gfa | 主要重叠群的GFA格式 |" >> "$report_file"
        echo "| Primary contigs (FASTA) | ${PREFIX}_PCTG.fasta.gz | 主要重叠群的FASTA格式（推荐使用） |" >> "$report_file"
        echo "| Alternate contigs | ${PREFIX}_hifiasm_run001.bp.a_ctg.gfa | 替代重叠群 |" >> "$report_file"
        echo "| Processed unitigs | ${PREFIX}_hifiasm_run001.bp.p_utg.gfa | 处理后的unitigs |" >> "$report_file"
        echo "| Raw unitigs | ${PREFIX}_hifiasm_run001.bp.r_utg.gfa | 原始unitigs |" >> "$report_file"
        echo "| Contig links | ${PREFIX}_PCTG.links | 重叠群连接信息 |" >> "$report_file"
        echo "" >> "$report_file"
        
        # 检查各种GFA文件并显示基本信息
        echo "#### 各类型组装文件统计" >> "$report_file"
        echo "" >> "$report_file"
        
        local gfa_files=(
            "${PREFIX}_hifiasm_run001.bp.p_ctg.gfa:Primary contigs"
            "${PREFIX}_hifiasm_run001.bp.a_ctg.gfa:Alternate contigs"  
            "${PREFIX}_hifiasm_run001.bp.p_utg.gfa:Processed unitigs"
            "${PREFIX}_hifiasm_run001.bp.r_utg.gfa:Raw unitigs"
        )
        
        echo "| 文件类型 | 序列数量 | 总长度 (bp) |" >> "$report_file"
        echo "|---------|---------|-------------|" >> "$report_file"
        
        for gfa_info in "${gfa_files[@]}"; do
            local gfa_file="${gfa_info%%:*}"
            local gfa_desc="${gfa_info##*:}"
            local full_path="${PREFIX}_results/$gfa_file"
            
            if [[ -f "$full_path" ]]; then
                local seq_count=$(grep -c "^S" "$full_path" 2>/dev/null || echo "0")
                local total_length=$(awk '/^S/ {sum+=length($3)} END {print sum+0}' "$full_path" 2>/dev/null || echo "0")
                echo "| $gfa_desc | $seq_count | $total_length |" >> "$report_file"
            else
                echo "| $gfa_desc | 文件不存在 | - |" >> "$report_file"
            fi
        done
    fi
    
    # 最后添加注意事项
    echo "" >> "$report_file"
    echo "**注意**: 流水线运行完成后，只保留 \`${PREFIX}_results/\` 目录，所有中间文件和.bin文件已被清理。" >> "$report_file"
    
    log_info "Markdown报告已保存到: $report_file"
}

# 整理hifiasm文件到results目录
organize_hifiasm_files() {
    log_step "整理Hifiasm文件到results目录"
    
    # 将有用的hifiasm文件移动到results目录（不包括.bin文件）
    log_info "移动hifiasm文件到results目录..."
    if [[ -d "${PREFIX}_temp/assembly_files" ]]; then
        # 移动所有非.bin文件
        find "${PREFIX}_temp/assembly_files" -type f \( -name "*.gfa" -o -name "*.fa" -o -name "*.paf" -o -name "*.log" -o -name "*.bed" \) | \
        while read file; do
            if [[ -f "$file" ]]; then
                mv "$file" "${PREFIX}_results/"
                log_info "移动文件: $(basename "$file")"
            fi
        done
        
        # 删除.bin文件
        log_info "删除hifiasm .bin文件..."
        find "${PREFIX}_temp/assembly_files" -name "*.bin" -type f -delete 2>/dev/null
        log_info ".bin文件已删除"
    fi
    
    log_info "Hifiasm文件整理完成"
}

# 清理临时文件和logs目录
cleanup() {
    log_step "清理临时文件和logs目录"
    
    # 备份k-mer直方图文件到results目录（用于报告显示）
    if [[ -f "${PREFIX}_temp/kmc_files/${PREFIX}_${KMER_SIZE}kmer.histo" ]]; then
        mkdir -p "${PREFIX}_results/kmc_files"
        cp "${PREFIX}_temp/kmc_files/${PREFIX}_${KMER_SIZE}kmer.histo" "${PREFIX}_results/kmc_files/"
        log_info "K-mer直方图文件已备份到results目录"
    fi
    
    # 备份软件版本信息到results目录
    if [[ -f "${PREFIX}_temp/software_versions.txt" ]]; then
        cp "${PREFIX}_temp/software_versions.txt" "${PREFIX}_results/"
        log_info "软件版本信息已备份到results目录"
    fi
    
    # 删除整个临时目录
    if [[ -d "${PREFIX}_temp" ]]; then
        rm -rf "${PREFIX}_temp"
        log_info "临时目录已删除: ${PREFIX}_temp"
    fi
    
    # 删除logs目录
    if [[ -d "${PREFIX}_logs" ]]; then
        rm -rf "${PREFIX}_logs"
        log_info "日志目录已删除: ${PREFIX}_logs"
    fi
    
    log_info "清理完成，只保留 ${PREFIX}_results 目录"
}

# 主函数
main() {
    log_info "开始基因组组装流水线"
    log_info "输入文件: $INPUT_BAM"
    log_info "输出前缀: $PREFIX"
    
    # 检查工具
    check_tools
    
    # 创建目录
    create_dirs
    
    # 记录参数
    log_parameters
    
    # 执行流水线步骤
    step1_data_extraction
    step2_genome_survey
    step3_genome_assembly
    
    # 整理hifiasm文件到results目录
    organize_hifiasm_files
    
    # 在清理之前先生成报告（需要access某些临时文件）
    generate_report
    
    # 清理临时文件
    cleanup
    
    log_info "基因组组装流水线完成！"
    log_info "所有结果保存在: ${PREFIX}_results/"
    log_info "最终报告: ${PREFIX}_results/${PREFIX}_assembly_report.md"
    log_info "其他临时文件已清理，只保留results目录"
}

# 运行主函数
main "$@"

